Setup: Libraries and utilities

import single cell data

data <- read_h5ad(here::here("data", "preprocessedData",  "GSE193816", "GSE193816_all_cells_data.h5ad"))
data <- CreateSeuratObject(counts = t(as.matrix(data$X)), 
                           meta.data = data$obs)

## mnp
mnp <- read_h5ad(here::here("data","preprocessedData",  "GSE193816", "GSE193816_mnp_data.h5ad"))
mnp <- CreateSeuratObject(counts = t(as.matrix(mnp$X)), 
                           meta.data = mnp$obs)

## t cells
tcell <- read_h5ad(here::here("data","preprocessedData",  "GSE193816", "GSE193816_t_cell_data.h5ad"))
tcell <- CreateSeuratObject(counts = t(as.matrix(tcell$X)), 
                          meta.data = tcell$obs)

quality metrics

All cells

VlnPlot(data, features = c("nFeature_RNA", "nCount_RNA", 
                           "percent_mito"), ncol = 5)

MNP cells

VlnPlot(mnp, features = c("nFeature_RNA", "nCount_RNA", 
                           "percent_mito"), ncol = 5)

T cells

VlnPlot(tcell, features = c("nFeature_RNA", "nCount_RNA", 
                           "percent_mito"), ncol = 5)

Dimension reduction plots

All cells

data <- NormalizeData(data,
                      normalization.method = "LogNormalize") %>% 
  FindVariableFeatures(selection.method = "vst", 
                       nfeatures = 2000) %>% 
  ScaleData(features = rownames(data)) %>% 
  RunPCA() %>% 
  RunUMAP(dims = 1:50, reduction.name = 'umap.rna', 
          reduction.key = 'rnaUMAP_')
Idents(data) <- data@meta.data$cluster

allcells <- DimPlot(data, reduction = "umap.rna", label = TRUE, 
                   label.size = 4, repel = TRUE, group.by = "cluster") + 
  ggtitle("RNA - all cells") +
  xlab("UMAP 1") +
  ylab("UMAP 2")
ggsave(here::here("results", "allcells_scrnaseq.png"), allcells, height = 4, width = 6)
allcells

save markers of all cells
res <- Seurat::FindAllMarkers(data, only.pos = TRUE)
cell_markers <- res %>% 
  group_by(cluster) %>% 
  dplyr::slice(1:10)
cellmarkers <- split(cell_markers$gene, cell_markers$cluster)
saveRDS(cellmarkers, here::here("results", "cellmarkers_GSE193816.rds"))

cellmarkers
## $`B cells`
##  [1] "MS4A1"     "BANK1"     "CD79A"     "HLA-DQA1"  "CD79B"     "TNFRSF13C"
##  [7] "HLA-DQB1"  "IGKC"      "POU2F2"    "HLA-DPB1" 
## 
## $`CD4 T cells`
##  [1] "LTB"    "CD3D"   "IL32"   "TRAC"   "IL7R"   "CD2"    "SPOCK2" "CD3E"  
##  [9] "CD3G"   "IKZF1" 
## 
## $`CD8 T cells`
##  [1] "CCL5" "CD3D" "CD2"  "CD3G" "CD8A" "IL32" "TRAC" "NKG7" "CD3E" "GZMA"
## 
## $`Epithelial cells`
##  [1] "KRT19"   "FXYD3"   "AGR2"    "WFDC2"   "ALDH1A1" "PERP"    "SLPI"   
##  [8] "KRT18"   "CXCL17"  "TACSTD2"
## 
## $MNP
##  [1] "TYROBP"   "FCER1G"   "AIF1"     "HLA-DQB1" "HLA-DQA1" "LST1"    
##  [7] "HLA-DPB1" "CCDC88A"  "MS4A6A"   "IGSF6"   
## 
## $`Mast cells`
##  [1] "GATA2"   "CPA3"    "SAMSN1"  "AKAP12"  "HDC"     "IL1RL1"  "MS4A2"  
##  [8] "CSF2RB"  "SLC18A2" "SLC45A3"
## 
## $`NK cells`
##  [1] "TRDC"      "XCL1"      "TNFRSF18"  "GNLY"      "XCL2"      "TXK"      
##  [7] "KLRF1"     "LINC00996" "SPINK2"    "SH2D1B"

MNP cells

mnp <- NormalizeData(mnp,
                      normalization.method = "LogNormalize") %>% 
  FindVariableFeatures(selection.method = "vst", 
                       nfeatures = 2000) %>% 
  ScaleData(features = rownames(mnp)) %>% 
  RunPCA() %>% 
  RunUMAP(dims = 1:50, reduction.name = 'umap.rna', 
          reduction.key = 'rnaUMAP_')

mnpcells <- DimPlot(mnp, reduction = "umap.rna", label = TRUE, 
                   label.size = 2.5, repel = TRUE, group.by = "annotation") + 
  ggtitle("RNA - MNP cells") +
  xlab("UMAP 1") +
  ylab("UMAP 2")
ggsave(here::here("results", "mnpcells_scrnaseq.png"), mnpcells, height = 4, width = 6)
mnpcells

T cells

tcell <- NormalizeData(tcell,
                      normalization.method = "LogNormalize") %>% 
  FindVariableFeatures(selection.method = "vst", 
                       nfeatures = 2000) %>% 
  ScaleData(features = rownames(tcell)) %>% 
  RunPCA() %>% 
  RunUMAP(dims = 1:50, reduction.name = 'umap.rna', 
          reduction.key = 'rnaUMAP_')

tcells <- DimPlot(tcell, reduction = "umap.rna", label = TRUE, 
                   label.size = 2.5, repel = TRUE, group.by = "annotation") + 
  ggtitle("RNA - T cells") +
  xlab("UMAP 1") +
  ylab("UMAP 2")
ggsave(here::here("results", "tcells_scrnaseq.png"), tcells, height = 4, width = 6)
tcells

Cell-type proportions

All cells

data@meta.data %>% 
  dplyr::select(id, group, condition, cluster) %>%
  dplyr::group_by(id, group, condition, cluster) %>% 
  mutate(n = n()) %>% 
  mutate(condition = factor(as.character(condition), 
                            levels = c("Bln", "Dil", "Ag"))) %>% 
  ggplot(aes(x = condition, y = n, color = group)) +
  geom_boxplot() +
  facet_wrap(~cluster, scales = "free") +
  theme_bw()

MNP cells

mnp@meta.data %>% 
  dplyr::select(id, group, condition, annotation) %>%
  dplyr::group_by(id, group, condition, annotation) %>% 
  mutate(n = n()) %>% 
  mutate(condition = factor(as.character(condition), 
                            levels = c("Bln", "Dil", "Ag"))) %>% 
  ggplot(aes(y = annotation, x = n, color = condition)) +
  geom_boxplot() +
  facet_grid(~group) +
  theme_bw()

T cells

tcell@meta.data %>% 
  dplyr::select(id, group, condition, annotation) %>%
  dplyr::group_by(id, group, condition, annotation) %>% 
  mutate(n = n()) %>% 
  mutate(condition = factor(as.character(condition), 
                            levels = c("Bln", "Dil", "Ag"))) %>% 
  ggplot(aes(y = annotation, x = n, color = condition)) +
  geom_boxplot() +
  facet_grid(~group) +
  theme_bw()

Differential expression analysis

All cells

sct_counts <- GetAssayData(data, assay = "RNA", layer = "counts")
metadata_aa <- data@meta.data %>% 
  dplyr::filter(group == "AA",
                condition != "Dil")
metadata_aa$group_id <- droplevels(metadata_aa$condition)
metadata_aa$sample_id <- droplevels(metadata_aa$Channel)
all_aa = sc_diffexp(sct_counts, metadata_aa, cell_ann_col = "cluster")
##   |                                                                              |                                                                      |   0%  |                                                                              |==========                                                            |  14%  |                                                                              |====================                                                  |  29%  |                                                                              |==============================                                        |  43%  |                                                                              |========================================                              |  57%  |                                                                              |==================================================                    |  71%  |                                                                              |============================================================          |  86%  |                                                                              |======================================================================| 100%
all_aa$toptables$condition <- "AA"
all_aa$toptables$cc <- "All cells"
metadata_ac <- data@meta.data %>% 
  dplyr::filter(group == "AC",
                condition != "Dil")
metadata_ac$group_id <- droplevels(metadata_ac$condition)
metadata_ac$sample_id <- droplevels(metadata_ac$Channel)
all_ac = sc_diffexp(sct_counts, metadata_ac, cell_ann_col = "cluster")
##   |                                                                              |                                                                      |   0%  |                                                                              |==========                                                            |  14%  |                                                                              |====================                                                  |  29%  |                                                                              |==============================                                        |  43%  |                                                                              |========================================                              |  57%  |                                                                              |==================================================                    |  71%  |                                                                              |============================================================          |  86%  |                                                                              |======================================================================| 100%
all_ac$toptables$condition <- "AC"
all_ac$toptables$cc <- "All cells"

saveRDS(all_aa$toptables,
        here::here("results", "sc_exacer_aa_toptable.rds"))

MNP cells

sct_counts <- GetAssayData(mnp, assay = "RNA", layer = "counts")
metadata_aa <- mnp@meta.data %>% 
  dplyr::filter(group == "AA",
                condition != "Dil")
metadata_aa$group_id <- droplevels(metadata_aa$condition)
metadata_aa$sample_id <- droplevels(metadata_aa$Channel)
mnp_aa = sc_diffexp(sct_counts, metadata_aa, cell_ann_col = "annotation")
##   |                                                                              |                                                                      |   0%  |                                                                              |=====                                                                 |   7%  |                                                                              |==========                                                            |  14%  |                                                                              |===============                                                       |  21%  |                                                                              |====================                                                  |  29%  |                                                                              |=========================                                             |  36%  |                                                                              |==============================                                        |  43%  |                                                                              |===================================                                   |  50%  |                                                                              |========================================                              |  57%  |                                                                              |=============================================                         |  64%  |                                                                              |==================================================                    |  71%  |                                                                              |=======================================================               |  79%  |                                                                              |============================================================          |  86%  |                                                                              |=================================================================     |  93%  |                                                                              |======================================================================| 100%
mnp_aa$toptables$condition <- "AA"
mnp_aa$toptables$cc <- "MNP cells"
metadata_ac <- mnp@meta.data %>% 
  dplyr::filter(group == "AC",
                condition != "Dil")
metadata_ac$group_id <- droplevels(metadata_ac$condition)
metadata_ac$sample_id <- droplevels(metadata_ac$Channel)
mnp_ac = sc_diffexp(sct_counts, metadata_ac, cell_ann_col = "annotation")
##   |                                                                              |                                                                      |   0%  |                                                                              |=====                                                                 |   7%  |                                                                              |==========                                                            |  14%  |                                                                              |===============                                                       |  21%  |                                                                              |====================                                                  |  29%  |                                                                              |=========================                                             |  36%  |                                                                              |==============================                                        |  43%  |                                                                              |===================================                                   |  50%  |                                                                              |========================================                              |  57%  |                                                                              |=============================================                         |  64%  |                                                                              |==================================================                    |  71%  |                                                                              |=======================================================               |  79%  |                                                                              |============================================================          |  86%  |                                                                              |=================================================================     |  93%  |                                                                              |======================================================================| 100%
mnp_ac$toptables$condition <- "AC"
mnp_ac$toptables$cc <- "MNP cells"

T cells

sct_counts <- GetAssayData(tcell, assay = "RNA", layer = "counts")
metadata_aa <- tcell@meta.data %>% 
  dplyr::filter(group == "AA",
                condition != "Dil")
metadata_aa$group_id <- droplevels(metadata_aa$condition)
metadata_aa$sample_id <- droplevels(metadata_aa$Channel)
tcell_aa = sc_diffexp(sct_counts, metadata_aa, cell_ann_col = "annotation")
##   |                                                                              |                                                                      |   0%  |                                                                              |=======                                                               |  10%  |                                                                              |==============                                                        |  20%  |                                                                              |=====================                                                 |  30%  |                                                                              |============================                                          |  40%  |                                                                              |===================================                                   |  50%  |                                                                              |==========================================                            |  60%  |                                                                              |=================================================                     |  70%  |                                                                              |========================================================              |  80%  |                                                                              |===============================================================       |  90%  |                                                                              |======================================================================| 100%
tcell_aa$toptables$condition <- "AA"
tcell_aa$toptables$cc <- "T cells"
metadata_ac <- tcell@meta.data %>% 
  dplyr::filter(group == "AC",
                condition != "Dil")
metadata_ac$group_id <- droplevels(metadata_ac$condition)
metadata_ac$sample_id <- droplevels(metadata_ac$Channel)
tcell_ac = sc_diffexp(sct_counts, metadata_ac, cell_ann_col = "annotation")
##   |                                                                              |                                                                      |   0%  |                                                                              |=======                                                               |  10%  |                                                                              |==============                                                        |  20%  |                                                                              |=====================                                                 |  30%  |                                                                              |============================                                          |  40%  |                                                                              |===================================                                   |  50%  |                                                                              |==========================================                            |  60%  |                                                                              |=================================================                     |  70%  |                                                                              |========================================================              |  80%  |                                                                              |===============================================================       |  90%  |                                                                              |======================================================================| 100%
tcell_ac$toptables$condition <- "AC"
tcell_ac$toptables$cc <- "T cells"

All cells

fdr_cutoff <- 0.2
diffexp <- rbind(all_ac$toptables, all_aa$toptables) %>%
  dplyr::rename(cell = cluster_id,
                adj.P.Val = p_adj.loc)

logfc = diffexp %>% 
  dplyr::select(gene, cell, logFC, cc, condition) %>% 
  spread(gene, logFC, fill = 0)
logfc2 <- as.matrix(logfc[, -c(1:4)])
rownames(logfc2) <- logfc$cell

cctype <- logfc$cc
cctype[logfc$cc == "All cells"] <- logfc$cell[logfc$cc == "All cells"]
cctype[cctype == "CD4 T cells"] <- "T cells"
cctype[cctype == "CD8 T cells"] <- "T cells"
cctype[cctype == "NK cells"] <- "T cells"
cctype[cctype == "MNP"] <- "MNP cells"

fdr = diffexp %>% 
  dplyr::select(gene, cell, adj.P.Val, cc, condition) %>% 
  spread(gene, adj.P.Val, fill = 1)
fdr2 <- fdr[, -c(1:4)]

row_ha <- rowAnnotation(cells=cctype, col = list(
  cells = c("B cells"="#999999", "Epithelial cells"="#009E73", 
            "Mast cells"="#F0E442", "MNP cells"="#D55E00", "T cells"="#56B4E9")
))
logfc$condition[logfc$condition == "AA"] <- "Allergic Asthma"
logfc$condition[logfc$condition == "AC"] <- "Allergic Controls"

ht_allcells <- Heatmap(logfc2, right_annotation = row_ha,
        row_split = logfc$condition,
        cluster_rows = FALSE,
        border = TRUE, name="logFC",
    cell_fun = function(j, i, x, y, width, height, fill) {
        if(fdr2[i, j] < fdr_cutoff)
            grid.text("*", x, y, gp = gpar(fontsize = 20))
})
ht_allcells

ht_allcells <- grid.grabExpr(draw(ht_allcells, merge_legends = FALSE))
ggsave(here("results", "all_cells_heatmap.pdf"), 
       plot = ht_allcells, width = 18, height = 5)

## ncells
ncells_all <- data@meta.data %>% 
  dplyr::group_by(group, condition, cluster) %>% 
  summarise(Freq = n()) %>% 
  mutate(group_condition = paste(group, condition, sep="_")) %>% 
  dplyr::filter(condition != "Dil") %>% 
  ungroup() %>% 
  dplyr::select(-c("group", "condition")) %>% 
  spread(group_condition, Freq) %>% 
  dplyr::rename(annotation = "cluster")

p1 <- scale(ncells_all[,-1], center = FALSE,
      scale = colSums(ncells_all[,-1])) %>% 
  as.data.frame() %>% 
  mutate(annotation = as.character(ncells_all$annotation)) %>% 
  gather(cond, prop, -annotation) %>% 
  mutate(cond = factor(cond, c("AC_Bln", "AC_Ag", "AA_Bln", "AA_Ag"))) %>% 
  ggplot(aes(x = cond, y = prop, fill = annotation)) +
  geom_bar(stat="identity", color = "black") +
  xlab("condition") +
  ylab("Cell proportions")
ggsave(here::here("results", "cell_prop_allcells.png"), p1)

p1

## nsig genes
nsig_all <- diffexp %>% 
  mutate(fc = ifelse(logFC > 0, "up", "down")) %>% 
  dplyr::group_by(condition, cell, fc) %>% 
  dplyr::filter(adj.P.Val < fdr_cutoff) %>% 
  summarise(n = n()) %>% 
  mutate(cond_fc = paste(condition,fc, sep="_")) %>% 
  ungroup() %>% 
  dplyr::select(-c(condition, fc)) %>% 
  spread(cond_fc, n, fill= 0)

knitr::kable(nsig_all)
cell AA_down AA_up AC_down AC_up
B cells 6 26 1 0
CD4 T cells 10 44 2 2
CD8 T cells 10 56 2 1
Epithelial cells 3 26 0 0
Mast cells 2 6 2 6
MNP 7 36 3 3
NK cells 3 17 0 0

MNP cells

diffexp <- rbind(mnp_aa$toptables, mnp_ac$toptables) %>% 
  dplyr::rename(cell = cluster_id,
                adj.P.Val = p_adj.loc)

logfc = diffexp %>%  
  dplyr::select(gene, cell, logFC, cc, condition) %>% 
  spread(gene, logFC, fill = 0)
logfc2 <- as.matrix(logfc[, -c(1:4)])
rownames(logfc2) <- logfc$cell

cctype <- logfc$cell

fdr = diffexp %>%  
  dplyr::select(gene, cell, adj.P.Val, cc, condition) %>% 
  spread(gene, adj.P.Val, fill = 1)
fdr2 <- fdr[, -c(1:4)]

row_ha <- rowAnnotation(cells=cctype)
logfc$condition[logfc$condition == "AA"] <- "Allergic Asthma"
logfc$condition[logfc$condition == "AC"] <- "Allergic Controls"


mnp_cells <- Heatmap(logfc2, right_annotation = row_ha,
        row_split = logfc$condition,
        border = TRUE, name="logFC",
    cell_fun = function(j, i, x, y, width, height, fill) {
        if(fdr2[i, j] < fdr_cutoff)
            grid.text("*", x, y, gp = gpar(fontsize = 20))
})
mnp_cells

mnp_cells <- grid.grabExpr(draw(mnp_cells, merge_legends = FALSE))
ggsave(here::here("results", "mnp_cells_heatmap.pdf"), 
       plot = mnp_cells, width = 18, height = 5)

## ncells
ncells_mnp <- mnp@meta.data %>% 
  dplyr::group_by(group, condition, annotation) %>% 
  summarise(Freq = n()) %>% 
  mutate(group_condition = paste(group, condition, sep="_")) %>% 
  dplyr::filter(condition != "Dil") %>% 
  ungroup() %>% 
  dplyr::select(-c("group", "condition")) %>% 
  spread(group_condition, Freq)

p2 <- scale(ncells_mnp[,-1], center = FALSE,
      scale = colSums(ncells_mnp[,-1])) %>% 
  as.data.frame() %>% 
  mutate(annotation = as.character(ncells_mnp$annotation)) %>% 
  gather(cond, prop, -annotation) %>% 
  mutate(cond = factor(cond, c("AC_Bln", "AC_Ag", "AA_Bln", "AA_Ag"))) %>% 
  ggplot(aes(x = cond, y = prop, fill = annotation)) +
  geom_bar(stat="identity", color = "black") +
  xlab("condition") +
  ylab("Cell proportions")
ggsave(here::here("results", "cell_prop_mnp.png"), p2, height = 5)

p2

## nsig genes
nsig_mnp <- diffexp %>% 
  mutate(fc = ifelse(logFC > 0, "up", "down")) %>% 
  dplyr::group_by(condition, cell, fc) %>% 
  dplyr::filter(adj.P.Val < fdr_cutoff) %>% 
  summarise(n = n()) %>% 
  mutate(cond_fc = paste(condition,fc, sep="_")) %>% 
  ungroup() %>% 
  dplyr::select(-c(condition, fc)) %>% 
  spread(cond_fc, n, fill= 0)
knitr::kable(nsig_mnp)
cell AA_down AA_up AC_down AC_up
DC2 (CD1C) 0 0 0 1
Mac1 (FABP4) 0 0 5 1
Mac2 (A2M) 2 2 0 0
MC2 (SPP1) 5 7 0 6
MC3 (AREG) 1 0 0 0
migDC (CCR7) 1 0 0 0
quiesMC 0 4 2 8

T cells

diffexp <- rbind(tcell_ac$toptables, tcell_aa$toptables) %>% 
  dplyr::rename(cell = cluster_id,
                adj.P.Val = p_adj.loc)

logfc = diffexp %>% 
  dplyr::select(gene, cell, logFC, cc, condition) %>% 
  spread(gene, logFC, fill = 0)
logfc2 <- as.matrix(logfc[, -c(1:4)])
rownames(logfc2) <- logfc$cell

cctype <- logfc$cell

fdr = diffexp %>% 
  dplyr::select(gene, cell, adj.P.Val, cc, condition) %>% 
  spread(gene, adj.P.Val, fill = 1)
fdr2 <- fdr[, -c(1:4)]
# fdr2[is.na(fdr2)] <- 1

row_ha <- rowAnnotation(cells=cctype)
logfc$condition[logfc$condition == "AA"] <- "Allergic Asthma"
logfc$condition[logfc$condition == "AC"] <- "Allergic Controls"

t_cells <- Heatmap(logfc2, right_annotation = row_ha, 
        row_split = logfc$condition,
        border = TRUE, name="logFC",
    cell_fun = function(j, i, x, y, width, height, fill) {
        if(fdr2[i, j] < fdr_cutoff)
            grid.text("*", x, y, gp = gpar(fontsize = 20))
})
t_cells

t_cells <- grid.grabExpr(draw(t_cells, merge_legends = FALSE))
ggsave(here::here("results", "t_cells_heatmap.pdf"), 
       plot = t_cells, width = 18, height = 5)

## ncells
ncells_tcells <- tcell@meta.data %>% 
  dplyr::group_by(group, condition, annotation) %>% 
  summarise(Freq = n()) %>% 
  mutate(group_condition = paste(group, condition, sep="_")) %>% 
  dplyr::filter(condition != "Dil") %>% 
  ungroup() %>% 
  dplyr::select(-c("group", "condition")) %>% 
  spread(group_condition, Freq)

p3 <- scale(ncells_tcells[,-1], center = FALSE,
      scale = colSums(ncells_tcells[,-1])) %>% 
  as.data.frame() %>% 
  mutate(annotation = as.character(ncells_tcells$annotation)) %>% 
  gather(cond, prop, -annotation) %>% 
  mutate(cond = factor(cond, c("AC_Bln", "AC_Ag", "AA_Bln", "AA_Ag"))) %>% 
  ggplot(aes(x = cond, y = prop, fill = annotation)) +
  geom_bar(stat="identity", color = "black") +
  xlab("condition") +
  ylab("Cell proportions")
ggsave(here::here("results", "cell_prop_tcells.png"), p3, height = 5)
p4 <- cowplot::plot_grid(p1, p2, p3, labels = c("A", "B", "C"), ncol = 1)
ggsave(here::here("results", "cell_prop_combined.png"), p4, height = 12)

p4

## nsig genes
nsig_tcells <- diffexp %>% 
  mutate(fc = ifelse(logFC > 0, "up", "down")) %>% 
  dplyr::group_by(condition, cell, fc) %>% 
  dplyr::filter(adj.P.Val < fdr_cutoff) %>% 
  summarise(n = n()) %>% 
  mutate(cond_fc = paste(condition,fc, sep="_")) %>% 
  ungroup() %>% 
  dplyr::select(-c(condition, fc)) %>% 
  spread(cond_fc, n, fill= 0)

knitr::kable(nsig_tcells)
cell AA_down AA_up AC_down AC_up
CD4 T (CD40LG) 8 33 0 0
CD4 Th17 (RORA) 5 37 2 2
CD4 Th2 (GATA3) 1 6 0 0
CD4 Treg (FOXP3) 0 0 1 0
CD8 T (CLIC3) 4 40 5 11
CD8 T (GZMK) 6 13 0 2
quiesCD8 T 4 9 0 0
Tgd (TRDC) 8 19 1 4

number of cells

ncells <- cbind(cellclass = c(as.character(ncells_all$annotation), rep("MNP", nrow(ncells_mnp)), rep("T cells", nrow(ncells_tcells))), 
                rbind(ncells_all, ncells_mnp, ncells_tcells))
write.csv(ncells[, c("cellclass", "annotation", "AC_Bln", "AC_Ag", "AA_Bln", "AA_Ag")], here::here("results", "ncells.csv"))

knitr::kable(ncells)
cellclass annotation AA_Ag AA_Bln AC_Ag AC_Bln
B cells B cells 582 434 184 173
CD4 T cells CD4 T cells 2598 531 579 1180
CD8 T cells CD8 T cells 1783 1993 1748 3583
Epithelial cells Epithelial cells 2387 5028 1616 7720
MNP MNP 3213 499 2112 946
Mast cells Mast cells 276 125 102 40
NK cells NK cells 189 60 93 94
MNP AS DC (AXL) 45 3 29 8
MNP Cycling (PCLAF) 3 16 8 18
MNP DC1 (CLEC9A) 9 2 8 7
MNP DC2 (CD1C) 557 24 183 43
MNP MC1 (CXCL10) 22 7 56 18
MNP MC2 (SPP1) 174 25 414 91
MNP MC3 (AREG) 188 34 373 76
MNP MC4 (CCR2) 895 14 326 20
MNP Mac1 (FABP4) 15 129 67 293
MNP Mac2 (A2M) 102 20 50 51
MNP migDC (CCR7) 100 21 77 30
MNP pDC (TCF4) 447 30 269 35
MNP quiesMC 541 148 210 159
MNP quiesMac 12 10 18 49
T cells CD4 T (CD40LG) 1570 330 369 873
T cells CD4 Th2 (GATA3) 199 33 33 23
T cells CD4 Th17 (RORA) 344 241 177 441
T cells CD4 ThIFNR (ISG15) 47 10 15 20
T cells CD4 Treg (FOXP3) 579 43 86 122
T cells CD8 T (CLIC3) 699 658 817 1403
T cells CD8 T (EGR2) 19 21 14 44
T cells CD8 T (GZMK) 154 82 160 306
T cells Tgd (TRDC) 360 305 245 295
T cells quiesCD8 T 377 545 388 956

number of significant gene-transcripts

nsig <- cbind(cellclass = c(as.character(nsig_all$cell), 
                       rep("MNP", nrow(nsig_mnp)), 
                       rep("T cells", nrow(nsig_tcells))), 
                rbind(nsig_all, nsig_mnp, nsig_tcells))
write.csv(nsig[, c("cellclass", "cell", "AC_up", "AC_down","AA_up", "AA_down")], here::here("results", "nsig.csv"))

knitr::kable(nsig)
cellclass cell AA_down AA_up AC_down AC_up
B cells B cells 6 26 1 0
CD4 T cells CD4 T cells 10 44 2 2
CD8 T cells CD8 T cells 10 56 2 1
Epithelial cells Epithelial cells 3 26 0 0
Mast cells Mast cells 2 6 2 6
MNP MNP 7 36 3 3
NK cells NK cells 3 17 0 0
MNP DC2 (CD1C) 0 0 0 1
MNP Mac1 (FABP4) 0 0 5 1
MNP Mac2 (A2M) 2 2 0 0
MNP MC2 (SPP1) 5 7 0 6
MNP MC3 (AREG) 1 0 0 0
MNP migDC (CCR7) 1 0 0 0
MNP quiesMC 0 4 2 8
T cells CD4 T (CD40LG) 8 33 0 0
T cells CD4 Th17 (RORA) 5 37 2 2
T cells CD4 Th2 (GATA3) 1 6 0 0
T cells CD4 Treg (FOXP3) 0 0 1 0
T cells CD8 T (CLIC3) 4 40 5 11
T cells CD8 T (GZMK) 6 13 0 2
T cells quiesCD8 T 4 9 0 0
T cells Tgd (TRDC) 8 19 1 4

sessionInfo

sessionInfo()
## R version 4.5.1 (2025-06-13)
## Platform: x86_64-apple-darwin20
## Running under: macOS Ventura 13.6.7
## 
## Matrix products: default
## BLAS:   /Library/Frameworks/R.framework/Versions/4.5-x86_64/Resources/lib/libRblas.0.dylib 
## LAPACK: /Library/Frameworks/R.framework/Versions/4.5-x86_64/Resources/lib/libRlapack.dylib;  LAPACK version 3.12.1
## 
## locale:
## [1] en_US.UTF-8/en_US.UTF-8/en_US.UTF-8/C/en_US.UTF-8/en_US.UTF-8
## 
## time zone: America/Vancouver
## tzcode source: internal
## 
## attached base packages:
## [1] grid      stats4    stats     graphics  grDevices utils     datasets 
## [8] methods   base     
## 
## other attached packages:
##  [1] here_1.0.2                  circlize_0.4.16            
##  [3] ComplexHeatmap_2.26.0       limma_3.66.0               
##  [5] scater_1.38.0               scuttle_1.20.0             
##  [7] muscat_1.24.0               SingleCellExperiment_1.32.0
##  [9] SummarizedExperiment_1.40.0 Biobase_2.70.0             
## [11] GenomicRanges_1.62.0        Seqinfo_1.0.0              
## [13] IRanges_2.44.0              S4Vectors_0.48.0           
## [15] BiocGenerics_0.56.0         generics_0.1.4             
## [17] MatrixGenerics_1.22.0       matrixStats_1.5.0          
## [19] lubridate_1.9.4             forcats_1.0.1              
## [21] stringr_1.6.0               dplyr_1.1.4                
## [23] purrr_1.1.0                 readr_2.1.6                
## [25] tidyr_1.3.1                 tibble_3.3.0               
## [27] ggplot2_4.0.1               tidyverse_2.0.0            
## [29] Seurat_5.3.1                SeuratObject_5.2.0         
## [31] sp_2.2-0                    anndata_0.8.0              
## 
## loaded via a namespace (and not attached):
##   [1] bitops_1.0-9             spatstat.sparse_3.1-0    httr_1.4.7              
##   [4] RColorBrewer_1.1-3       doParallel_1.0.17        numDeriv_2016.8-1.1     
##   [7] backports_1.5.0          tools_4.5.1              sctransform_0.4.2       
##  [10] R6_2.6.1                 lazyeval_0.2.2           uwot_0.2.4              
##  [13] mgcv_1.9-3               GetoptLong_1.0.5         withr_3.0.2             
##  [16] prettyunits_1.2.0        gridExtra_2.3            fdrtool_1.2.18          
##  [19] progressr_0.18.0         textshaping_1.0.4        cli_3.6.5               
##  [22] Cairo_1.7-0              spatstat.explore_3.5-3   fastDummies_1.7.5       
##  [25] sandwich_3.1-1           labeling_0.4.3           slam_0.1-55             
##  [28] sass_0.4.10              mvtnorm_1.3-3            S7_0.2.0                
##  [31] spatstat.data_3.1-9      blme_1.0-6               ggridges_0.5.7          
##  [34] pbapply_1.7-4            systemfonts_1.3.1        parallelly_1.45.1       
##  [37] rstudioapi_0.17.1        shape_1.4.6.1            gtools_3.9.5            
##  [40] ica_1.0-3                spatstat.random_3.4-2    Matrix_1.7-3            
##  [43] ggbeeswarm_0.7.2         abind_1.4-8              lifecycle_1.0.4         
##  [46] yaml_2.3.12              edgeR_4.8.0              gplots_3.2.0            
##  [49] SparseArray_1.10.1       Rtsne_0.17               promises_1.5.0          
##  [52] crayon_1.5.3             miniUI_0.1.2             lattice_0.22-7          
##  [55] beachmat_2.26.0          cowplot_1.2.0            magick_2.9.0            
##  [58] pillar_1.11.1            knitr_1.50               rjson_0.2.23            
##  [61] boot_1.3-31              estimability_1.5.1       corpcor_1.6.10          
##  [64] future.apply_1.20.0      codetools_0.2-20         glue_1.8.0              
##  [67] spatstat.univar_3.1-4    data.table_1.17.8        vctrs_0.6.5             
##  [70] png_0.1-8                spam_2.11-1              Rdpack_2.6.4            
##  [73] gtable_0.3.6             assertthat_0.2.1         cachem_1.1.0            
##  [76] xfun_0.54                rbibutils_2.4            S4Arrays_1.10.0         
##  [79] mime_0.13                coda_0.19-4.1            reformulas_0.4.2        
##  [82] survival_3.8-3           iterators_1.0.14         statmod_1.5.1           
##  [85] fitdistrplus_1.2-4       ROCR_1.0-11              nlme_3.1-168            
##  [88] pbkrtest_0.5.5           EnvStats_3.1.0           progress_1.2.3          
##  [91] RcppAnnoy_0.0.22         rprojroot_2.1.1          bslib_0.9.0             
##  [94] TMB_1.9.18               irlba_2.3.5.1            vipor_0.4.7             
##  [97] KernSmooth_2.23-26       otel_0.2.0               colorspace_2.1-2        
## [100] ggrastr_1.0.2            DESeq2_1.50.0            tidyselect_1.2.1        
## [103] emmeans_2.0.1            compiler_4.5.1           BiocNeighbors_2.4.0     
## [106] DelayedArray_0.36.0      plotly_4.11.0            caTools_1.18.3          
## [109] scales_1.4.0             remaCor_0.0.20           lmtest_0.9-40           
## [112] digest_0.6.37            goftest_1.2-3            presto_1.0.0            
## [115] spatstat.utils_3.2-0     minqa_1.2.8              variancePartition_1.40.0
## [118] rmarkdown_2.30           aod_1.3.3                RhpcBLASctl_0.23-42     
## [121] XVector_0.50.0           htmltools_0.5.8.1        pkgconfig_2.0.3         
## [124] lme4_1.1-37              lpsymphony_1.38.0        fastmap_1.2.0           
## [127] rlang_1.1.6              GlobalOptions_0.1.2      htmlwidgets_1.6.4       
## [130] shiny_1.11.1             farver_2.1.2             jquerylib_0.1.4         
## [133] IHW_1.38.0               zoo_1.8-14               jsonlite_2.0.0          
## [136] BiocParallel_1.44.0      BiocSingular_1.26.0      magrittr_2.0.4          
## [139] dotCall64_1.2            patchwork_1.3.2          Rcpp_1.1.0              
## [142] viridis_0.6.5            reticulate_1.44.0        stringi_1.8.7           
## [145] MASS_7.3-65              plyr_1.8.9               parallel_4.5.1          
## [148] listenv_0.10.0           ggrepel_0.9.6            deldir_2.0-4            
## [151] splines_4.5.1            tensor_1.5.1             hms_1.1.4               
## [154] locfit_1.5-9.12          igraph_2.2.1             spatstat.geom_3.6-0     
## [157] RcppHNSW_0.6.0           reshape2_1.4.4           ScaledMatrix_1.18.0     
## [160] evaluate_1.0.5           nloptr_2.2.1             tzdb_0.5.0              
## [163] foreach_1.5.2            httpuv_1.6.16            RANN_2.6.2              
## [166] polyclip_1.10-7          future_1.67.0            clue_0.3-66             
## [169] scattermore_1.2          rsvd_1.0.5               broom_1.0.10            
## [172] xtable_1.8-4             fANCOVA_0.6-1            RSpectra_0.16-2         
## [175] later_1.4.4              ragg_1.5.0               viridisLite_0.4.2       
## [178] lmerTest_3.1-3           glmmTMB_1.1.13           beeswarm_0.4.0          
## [181] cluster_2.1.8.1          timechange_0.3.0         globals_0.18.0